use "$dataDir/dataset_birthoutcomes_pollution", clear

local controls "edad_madre total_consultas"
local vd "weight gweeks vlowb vpreterm"

foreach vd in `vd'{
	reghdfe  `vd' pm25_pred  `controls', abs(geoid year) vce(cluster geoid)
		estimates store `vd'_pm25_pred
	reghdfe  `vd' l1_pm25_pred  `controls', abs(geoid year) vce(cluster geoid)
		estimates store `vd'_l1_pm25_pred
	reghdfe  `vd' l2_pm25_pred  `controls', abs(geoid year) vce(cluster geoid)
		estimates store `vd'_l2_pm25_pred
}


use "$dataDir/dataset_birthoutcomes_vulnerable_mothers_pollution", clear

local controls "edad_madre total_consultas"
local vd "weight gweeks vlowb vpreterm"

foreach vd in `vd'{
	reghdfe  `vd' pm25_pred  `controls', abs(geoid year) vce(cluster geoid)
		estimates store `vd'_pm25_pred_vm
	reghdfe  `vd' l1_pm25_pred  `controls', abs(geoid year) vce(cluster geoid)
		estimates store `vd'_l1_pm25_pred_vm
	reghdfe  `vd' l2_pm25_pred  `controls', abs(geoid year) vce(cluster geoid)
		estimates store `vd'_l2_pm25_pred_vm

}

la var pm25_pred "Pollution{sub:t}"
la var l1_pm25_pred "Pollution{sub:t-1}"
la var l2_pm25_pred "Pollution{sub:t-2}"

coefplot (weight_pm25_pred weight_l1_pm25_pred weight_l2_pm25_pred)  (weight_pm25_pred_vm weight_l1_pm25_pred_vm weight_l2_pm25_pred_vm) , keep(pm25_pred l1_pm25_pred l2_pm25_pred) vertical ciopts(recast(rcap))  yline(0, lcolor(black)) label p1(label("All"))  p2(label("Seguro Popular")) title("Panel a) Birth weight (g)") legen(cols(2)) saving("$tempDir/coefplot_birthweight.gph", replace)

coefplot (gweeks_pm25_pred gweeks_l1_pm25_pred gweeks_l2_pm25_pred)  (gweeks_pm25_pred_vm gweeks_l1_pm25_pred_vm gweeks_l2_pm25_pred_vm) , keep(pm25_pred l1_pm25_pred l2_pm25_pred) vertical ciopts(recast(rcap))  yline(0, lcolor(black)) label p1(label("All"))  p2(label("Seguro Popular")) title("Panel b) Gestation length (weeks)") legen(cols(2)) saving("$tempDir/coefplot_gweeks.gph", replace)

coefplot (vlowb_pm25_pred vlowb_l1_pm25_pred vlowb_l2_pm25_pred)  (vlowb_pm25_pred_vm vlowb_l1_pm25_pred_vm vlowb_l2_pm25_pred_vm), keep(pm25_pred l1_pm25_pred l2_pm25_pred)  vertical ciopts(recast(rcap))  rescale(pm25_pred = 100 l1_pm25_pred = 100 l2_pm25_pred  = 100) yline(0, lcolor(black)) p1(label("All"))  p2(label("Seguro Popular")) title("Panel c) Very low birth weight incidence") legen(cols(2)) saving("$tempDir/coefplot_vlowb.gph", replace)

coefplot (vpreterm_pm25_pred vpreterm_l1_pm25_pred vlowb_l2_pm25_pred)  (vpreterm_pm25_pred_vm vpreterm_l1_pm25_pred_vm vpreterm_l2_pm25_pred_vm), keep(pm25_pred l1_pm25_pred l2_pm25_pred)  vertical ciopts(recast(rcap))  rescale(pm25_pred = 100 l1_pm25_pred = 100 l2_pm25_pred  = 100) yline(0, lcolor(black))  p1(label("All"))  p2(label("Seguro Popular")) title("Panel d) Very preterm birth incidence") legen(cols(2)) saving("$tempDir/coefplot_vpreterm.gph", replace)
	
	grc1leg "$tempDir/coefplot_birthweight.gph"  "$tempDir/coefplot_gweeks.gph" "$tempDir/coefplot_vlowb.gph" "$tempDir/coefplot_vpreterm.gph", cols(2)
		graph export "$mainDir/figures/Fig6_main_health.pdf", as(pdf) replace